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Abstract 

A continuous Trigonometrically-fitted Second Derivative Method (CTSDM) whose coefficients depend on the 
frequency and stepsize is constructed using trigonometric basis functions. A discrete Trigonometrically-fitted second 
derivative method (TSDM) is recovered from the CTSDM as a by-product and applied to solve initial value problems 
(IVPs) with oscillating solutions. We discuss the stability properties of the TSDM and present numerical experiments to 
demonstrate the efficiency of the method. 
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Introduction 

In this paper, we consider the subclass of first order 
differential equation 



/ = f(x,y), y(a) = yo, x e [a, b] , 



(i) 



with periodic or oscillating solutions where / : $i x 
m m wn fytyQ e Oscillatory IVPs frequently arise 
in areas such as classical mechanics, celestial mechan- 
ics, quantum mechanics, and biological sciences. Several 
numerical methods based on the use of polynomial basis 
functions have been developed for solving this class of 
important problems (see Lambert , Hairer et al. in (Hairer 
and Wanner 1996), Hairer 1982, and Sommeijer 1993). 
Other methods based on exponential fitting techniques 
which take advantage of the special properties of the solu- 
tion that may be known in advance have been proposed 
(see Simos 1998, Vanden Berghe et al. 2001a, Vanden 
Berghe et al. 2009, Vigo-Aguiar et al. 2003, Franco 2002, 
Fang et al. 2009, Nguyen et al. 2007, Ozawa 2005, Jator 
et al. 2012, and Ngwane et al. 2012b). In the spirit of 
2005, the motivation governing the exponentially-fitted 
methods is inherent to the fact that if the frequency or a 
reasonable estimate of it is known in advance, these meth- 
ods will be more advantageous than the polynomial based 
methods. 
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The aim of this paper is to construct a TSDM. This 
construction is done by initially developing a CTSDM 
which then provides a discrete method that is applied as 
a TSDM which takes the frequency of the solution as a 
priori knowledge. In particular, CTSDM consists of a sum 
of continuous functions while TSDM is a by-product of 
CTSDM. The coefficients of the TSDM are functions of 
the frequency and the stepsize, hence the solutions pro- 
vided by the proposed method are highly accurate if (1) 
has periodic solutions with known frequencies. We adopt 
the approach given in Jator et al. in (Ngwane and Jator 
2012a; Jator et al. 2012), where the TSDM is used to obtain 
the approximation y n +\ to the exact solution y(x n -\-i) on 
the interval [x n ,x n +i], h — x n +\ — x ni n = 0, . . . ,N — 1, 
on a partition [a, b], where a, b e R, h is the constant step- 
size, n is a grid index and N > 0 is the number of steps. 
We note that second derivative methods with polynomial 
basis functions were proposed to overcome the Dahlquist 
1956 barrier theorem whereby the conventional linear 
multistep method was modified by incorporating the sec- 
ond derivative term in the derivation process in order to 
increase the order of the method, while preserving good 
stability properties (see Enright 1974). 

This paper is organized as follows. In Section "Develop- 
ment of method", we obtain a trigonometric basis repre- 
sentation U (x) for the exact solution which is used to 
generate a TSDM for solving (1). The analysis and imple- 
mentation of the TSDM are discussed in Section "Error 
analysis and stability". Numerical examples are given in 
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Section "Numerical examples" to show the accuracy and 
efficiency of the TSDM. Finally, we give some concluding 
remarks in Section "Conclusion". 

Development of method 

In this section, our objective is to construct a CTSDM 
which produces a discrete method as a by-product. The 
method has the form 



Jn+l — Jn 



(2) 



where u — wh, Pj(u), Yj( u )>j = 0, 1, are coefficients that 
depend on the stepsize and frequency. We note that y n +j 
is the numerical approximation to the analytical solution 

y(x n +j), and 

I x n+j 

Jn+j = J KXn+j, yn+j), gn+j = ^ \y n+j 

with ; = 0, 1. In order to obtain equation (2) we proceed 
by seeking to approximate the exact solution y(x) on the 
interval [x n ,x n + h] by the interpolating function U(x) of 
the form 

U (x) = ao+ a\x + 0L2X 2 + sin(wx) + cos(w#), (3) 

where a^a^a^a^ and are coefficients that must be 
uniquely determined. We then impose that the interpolat- 
ing function in (3) coincides with the analytical solution at 
the point x n to obtain the equation 



U (x n ) = y n - 



(4) 



We also demand that the function (3) satisfies the differ- 
ential equation (1) at the points x n +j, j = 0, 1 to obtain the 
following set of three equations: 



U' { x n+j) =fn+j, U" (x n +j) = gn+j, j = 0, 1. 



(5) 



Equations (4) and (5) lead to a system of five equations 
which is solved by Cramers rule to obtain aj, j = 
0, 1, 2, 3, 4. Our continuous CTSDM is constructed by 
substituting the values of aj into equation (3). After some 
algebraic manipulation, the CTSDM is expressed in the 
form 

U{x) =y n + h(po(w,x)f n + Pi(w,x)f n +i) + h 2 (y 0 (w,x)g n 
+ Yi(w,x)g n +i), (6) 

where w is the frequency, /3o(w,x), /3i(w,x), Yq(w,x), 
and yi(w,x), are continuous coefficients. The continuous 
method (6) is used to generate the method of the form (2). 



Thus, evaluating (6) at x = x n +i and letting u = wh, we 
obtain the coefficients of (2) as follows: 

Po = \, 
Pi = \, 

yo = (-csc(|)(^cos(|)-2sin(|)))/(2^), 

Yi = (esc (|) (u cos (f) - 2 sin (§))) / (2u 2 ) . 

(7) 

Error analysis and stability 

Local truncation error 

We note that when u 0 the coefficients given by (7) are 
vulnerable to heavy cancellations and hence the following 
Taylor series expansion must be used (see Simos 1998). 



1 _1_ u' 
Y0= 12 + 72 



720 ^ 30240 ^ 1209600 ^ 47900160 ^ 1307674368000 



Yi = - 



12 720 30240 1209600 47900160 1307674368000 



(8) 



In fact, for practical computations when u is small, it 
is better to use the series expansion (8) (see Calvo et al. 
2009). 

Thus the Local Truncation Error (LTE) of (2) subject to 
(8) is obtained as 

LTE = y(Xn+l) ~yn+l = ^ (x n ) +/ 5) (*„)) + O (h 6 ) . 

(9) 

Remark L The method (2) specified by (8) is a fourth- 
order method and reduces to a one-step conventional sec- 
ond derivative method as u -> 0 (see Lambert 1973, 
p. 201). 

Stability 

Proposition L The TSDM (2) applied to the test 
equations y' = ky andy" — X 2 y yields 

y n +\ = M(q; u)y n , q = hk, u = wh, (10) 

with 

M(q;u) = (l+qPo(u)+q 2 y 0 (u))~ 1 (l - qPi(u)-q 2 yi(u)). 

(ID 

Proof We begin by applying (2) to the test equations 
y = Xy andy r = X 2 y which are expressed 2isf(x,y) = ky 
and g(x, y) = X 2 y respectively; letting q = hX and u = wh, 
we obtain a linear equation which is used to solve for y n+ \ 
with (11) as a consequence. □ 
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Remark 2. The rational function M(q; u) is called the 
stability function which determines the stability of the 
method. 

Definition 1. A region of stability is a region in theq — u 
plane, in which \M(q; u) \ < 1. 



Remark 3. The TSDM (12) is consistent as it has order 
p > 1 and zero-stable, hence it is convergent since zero- 
stability + consistency = convergence. 

Corollary 1. The TSDM (12) has M(q; u) specified by 



The TSDM method (2) specified by (7) is given by 

y n +i =y n + -ifn + ^ csc (^/ 2 ) (12) 

x (cosO/2) - 2sin(u/2))(-g n +g n +{). 

Definition 2. The method (12) is zero-stable provided 
the roots of the first characteristic polynomial have modu- 
lus less than or equal to unity and those of modulus unity 
are simple (see Lambert 1991). 



M(q; u) = 1 + 



x / 1 



q csc(f)(^cos(f)-2sin(f))^ \ 
2 2u 2 J 

q csc(f) (wcos(f)-2sin(f))# 2N 



2u 2 



Remark 4. In the q — u plane the TSDM (12) is stable for 
q < 0, and u e[—2n, 2n], since from above \M(q; u)\ < 1, 
q < 0. 



Definition 3. The method (12) is consistent if it has 
order p > 1 (see (Fatunla 1991)). 



Remark 5. Figure 1 is a plot of the stability region and 
Figure 2 shows the zeros and poles ofM(q; u). We note from 



Oh 

I 



-100 < q < 100 



-100 -50 0 

Figure 1 The shaded region represents the truncated region of absolute stability. 



50 
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Ngwane and Jator SpringerPlus 201 4, 3:304 
http://www.springerplus.eom/content/3/1/304 



Page 4 of 1 1 






Figure 2 M(q; u) has zeros(D) and no poles(+) in C , with u = it. 



Figure 2 that the stability region includes the entire leftside 
of the complex plane. 

Definition 4. The TSDM with the stability function (11) 
is said to be A-stable at u = uo, — 2tt < uq < 2n, if 
\M(q; u)\ <l,Vq e C~, (see Nguyen etal. 2007). 

Remark 6. We observe from definition 1, remarks 4, 
5, and Figure 2, that TSDM is A-stable. In particular, 
\M(q; iy)\ — 1 Vy G R, and by the maximum principle, 



Table 1 Results with co = 1, e = 0.005, for Example 1 





TSDM 




FESDIRK4(3) 




ESDIRK4(3) 


N 


\Error\ 




N 


\Error\ 




N 


| Error | 


150 


1.203 x 10- 


-2 


170 


2.866 x 1 0" 


-1 


277 


2.153 x 10° 


200 


5.694 x 10" 


-3 


225 


7.846 x 1 0" 


-3 


496 


1.494 x 10~ 1 


300 


3.143 x 10- 


-4 


381 


1.399 x 10" 


-3 


884 


9.359 x 10- 3 


600 


1.259 x 10" 


-6 


680 


1 .690 x 1 0" 


-4 


1573 


6.200 x 10~ 4 


800 


1 .264 x 1 0" 


-7 


1207 


1 .846 x 1 0" 


-5 


2796 


4.416 x 10~ 5 


1600 


4.947 x 1 0" 


10 


2144 


1.938 x 10" 


-6 


4970 


3.412 x 10~ 6 


2400 


1.931 x 10- 


11 


3806 


1 .993 x 1 0" 


-7 


8833 


2.848 x 10~ 7 


3200 


1 .944 x 1 0" 


12 


6762 


2.021 x 10" 


-8 


15706 


2.530 x 10~ 8 



the method will be A-stable if\M(q; u) \ has no poles in the 
left plane (see E. Hairer et al. 1996, p. 43, 53). Moreover, the 
real part of the zeros of \M(q; u)\ must be negative, while 
the real part of the poles of\M(q; u) \ must be positive. 

Implementation 

In the spirit of Ngwane et al in (Ngwane and Jator 2012a; 
2012b), the TSDM (12) is implemented to solve (1) with- 
out requiring starting values and predictors. For instance, 
if we let n = 0 in (12), then y\ is obtained on the sub- 
interval [#o>#iL as Jo is known from the IVR Similarly, if 
n = 1, thenj2 is obtained on the sub-interval [#i,#2L as 
yi is known from the previous computation, and so on; 
until we reach the final sub-interval [xm-i>xn]- We note 



Table 2 Results with 


(O = 


1.01, for Example 2 




TSDM 




Simos 


Ixaru et al. 


N \Error\ 


N 


\Error\ 


N 


\Error\ 


150 3.3 x10- 3 


300 


1.7 x 10~ 3 


300 


1.1 x 10~ 3 


300 6.4 x10~ 5 


600 


1.9 x 10~ 4 


600 


5.4 x 10~ 5 


600 5.1 x 10~ 6 


1200 


1.4 x 10~ 5 


1200 


1.9 x 10~ 6 


2000 1.0 x10- 7 


2400 


8.7 x 10~ 7 


2400 


6.2 x 10- 8 
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Table 3 Results with co = 1 0, for Example 3 

TSDM Simos1998 



N 


\Error\ 


NFEs 


\Error\ 


NFEs 


1000 


1.7 x 10- 3 


4004 


1.4 x 10- 1 


8000 


2000 


2.5 x 10- 4 


8004 


3.5 x 10~ 2 


16000 


4000 


2.7 x 10- 5 


16004 


1.1 x 1CT 3 


32000 


8000 


1.6 x 1CT 6 


32004 


8.4 x 1CT 5 


64000 


16000 


1.0 x 10 -7 


64004 


5.5 x 10~ 6 


1 28000 


32000 


6.3 x 10- 9 


1 28004 


3.5 x 10- 7 


256000 



that for linear problems, we solve (1) directly using the 
feature solve[ ] in Matlab, while nonlinear problems use 
the Newton's method in Matlab enhanced by the feature 
fsolve[ ]. 

Numerical examples 

In this section, we give numerical examples to illustrate 
the accuracy (small errors) and efficiency (fewer number 
of function evaluations (NFEs)) of the TSDM. We find the 
approximate solution on the partition n^, where tt^ : a = 

Xo < X\ < X2 < ... < X n < X n +i < ... < Xn = 

b, and we give the errors at the endpoints calculated as 
Error =jjv — y(xx). We note that the method requires only 
two function evaluations per step and in general requires 
(2N + 2) NFEs on the entire interval. All computations 
were carried out using a written code in Matlab. 

Example 1. Consider the given two-body problem which 
was solved by Ozawa 200S. 

ji(0) = 1 - e, 34(0) = 0, y 2 (0) = 0, / 2 (0) 
[T+~e 

= J- , x g[0,50tt], 

V 1 - e 

where e, 0 < e < 1 is an eccentricity. The exact solution of 
this problem is 

Exact : y\(x) = cos(/<) — e, = \[\ — e 2 sin(/<), 

where k is the solution of the Kepler's equation k = x + 
e sin(/<). We choose co = 1. 



Table 4 Results with oo = 1, for Example 4 



TSDM Nguyen etal. 2007 



N 


\Error\ 


NFEs 


N 


\Error\ 


NFEs 


10 


1.3 x 1(T 15 


88 


73 


3.3 x 1CT 12 


327 


43 


8.4 x 10- 14 


368 


142 


0.9 x 10- 11 


707 


80 


7.1 x TO" 15 


648 


170 


3.7 x TO" 12 


811 



Table 5 Results with co = 314.16, for Example 5 on [0,100] 





TSDM 






CHEBY24 




N 


\Error\ 


NFEs 


N 


\Error\ 


NFEs 


9 


5.9 x 10- 14 


40 


9 


1.84 x 10- 11 


450 


20 


4.0 x TO" 15 


84 









Table 1 contains the results obtained using the TSDM. 
These results are compared with the explicit singly diago- 
nally implicit Runge-Kutta (ESDIRK) and the functionally 
fitted ESDIRK (FESDIRK) methods given in Ozawa 2005. 
In terms of accuracy, Table 1 clearly shows that TSDM 
performs better than those in Ozawa 2005. 

Example 2. We consider the nonlinear Duffing equation 
which was also solved Ixaru et al 2004. 

/' + y + y 3 =B cos(Qx), y(0) = C 0 , / (0) =0, xe[0, 300] . 
The analytic solution is given by 
Exact : y(x) = C\ cos(^) + C2 cos(3£2#) + C3 cos(5£2#) 
+ C4 cos (7 Six), 

where Q = 1.01, B = 0.002, C 0 = 0.200426728069, C\ = 
0.200179477536, C 2 = 0.246946143 x 10~ 3 , C 3 = 
0.304016 x 10~ 6 , C 4 = 0.374 x 10~ 9 . We choose co = 1.01 
and for more on frequency choice see Ramos et al. 2010. 

We compare the end-point global errors for TSDM with 
the fourth order methods in Ixaru et al. 2004. We see from 
Table 2 that the results produced by TSDM are better than 
Simos' method used in (Ixaru and Berghe 2004), as TSDM 
produces better error magnitude while using less number 
of steps and fewer number of function evaluations. TSDM 
is very competitive to the method used by Ixaru et al. 
2004. 

Example 3. We consider the following inhomogeneous 
IVP by Simos 1998. 

y " = -IOO3/ + 99sin(*), y(0) = 1, /(0) 
= 11, x e [0,1000] 
where the analytic solution is given by 

Exact : y(x) = cos(10^) + sin(10^) + sin(^). 

The exponentially-fitted method in Simos 1998 is fourth 
order and hence comparable to our method, TSDM. We 



Table 6 Results with co = 314.16, for Example 5 on [0, 1] 





TSDM 






CHEBY1 




N 


\Error\ 


NFEs 


N 


\Error\ 


NFEs 



1 1.29 x 1CT 21 8 1 1 x 1CT 16 8 
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Table 7 Results with oo = 1, for Example 6 with = -3 



TSDM with(/? = -3) Nguyen et al. 2007 with (/? = -3) 



N 


\Error\ 


NFEs 


N 


\Error\ 


NFEs 


6 


8.9 x 10~ 6 


28 


10 


5.4 x 10~ 6 


47 


10 


9.0 x 10~ 7 


44 


19 


8.3 x 10~ 8 


88 


27 


1.3 x 10~ 8 


112 


23 


4.5 x 10~ 4 


113 


40 


2.7 x 10~ 9 


164 









see from Table 3 that TSDM is more efficient than the 
method in Simos 1998. We also compare the computa- 
tional efficiency of the two methods by considering the 
NFEs over N integration steps for each method. Our 
method, TSDM, requires only 2N+2 function evaluations 
in N steps compared to 4N function evaluations in N steps 
for the method in Simos 1998. Hence for this example, 
TSDM performs better. 

Example 4. Linear Kramarz problem We consider the 
following second- or der IVP, (see Nguyen etal. 2007 [p. 204]) 

yf(0) = (o)' °^^ loa 

Exact : y(t) = (2 cos(£)> — cos(£)) r . 

We use this example to show the efficiency of TSDM on 
linear systems. Nguyen et al. 2007 used the "trigonometric 
implicit Runge-Kutta", TIRK3, method to solve the above 
linear Kramarz problem. Clearly, TSDM performs better 
as seen in Table 4. 

Example 5. We consider the IVP (see Vigo-Aguiar et al. 
2003) 

y" + K 2 y = K 2 x, y(0) = 10" 5 , / (0) 

= 1 - 7C10" 5 cot(/C), x e [0, 100] 

where K = 314.16, and we choose co = 314.16. The 
analytic solution is given by 

Exact : y(x) = x + 10 _5 (cos(/<#) — cot(7<T) sm(Kx)). 
Table 8 Results with oo = 1 , for Example 6 with j$ = -1000 



TSDM with (/? = - 1 000) Nguyen et al. 2007 with (/? = - 1 000) 



N 


\Error\ 


NFEs 


N 


\Error\ 


NFEs 


6 


8.9 x 10~ 6 


28 


13 


1.0 x 10~ 6 


61 


16 


1.2 x 10- 7 


68 


16 


1.6 x 10- 7 


76 


24 


2.3 x TO" 8 


100 


21 


7.0 x TO" 8 


98 



Table 9 Results, with predictor-corrector (PreCor) and 
co = 1 .01/ for Example 2 



TSDM 




PreCor 




Simos 




Ixaru et al. 




N \Error\ 


CPU 


\Error\ 


CPU 


N 


\Error\ 


N 


\Error\ 


150 3.3(-3) 


1.6 


1.7(-2) 


0.76 


300 


1.7(-3) 


300 


1.1 (-3) 


300 6.4(-5) 


2.3 


4.0(-4) 


1.6 


600 


1.9(-4) 


600 


5.4(-5) 


600 5.1 (-6) 


5.5 


1.2(-4) 


2.9 


1200 


1.4(-5) 


1200 


1.9(-6) 


2000 1.0(-7) 


18.2 


2.0(-5) 


10 


2400 


87(-7) 


2400 


6.2(-8) 



This problem demonstrates the performance of TSDM 
on a well-known oscillatory problem. We compare the 
results from TSDM with the Dissipative Chebyshev 
exponential-fitted methods, CHEBY24 and CHEBY1 used 
in Vigo-Aguiar et al. 2003. We see that TSDM uses fewer 
number of function evaluations with better accuracy than 
CHEBY24 that is designed to use fewer number of steps. 
Integrating in the interval [0, 1] with a stepsize equal to 
the total length of the interval, we obtain an error of 
order 10 -21 . Hence TSDM is a more efficient integra- 
tor. We note that compared with the methods CHEBY24 
and CHEBY1 which use stepsizes considerably larger than 
those used in multistep methods, TSDM is very competi- 
tive and superior to both CHEBY24 and CHEBY1. 

Example 6. A nearly sinusoidal problem 
We consider the following IVP on the range 0 < t < 10, 
(see Nguyen et al 2007, p. 205) 

y[ = -2yi +y 2 + sin(f), yi(0) = 2 

/ 2 = -(^+2)ji + (^ + l)j 2 +sin(0-cos(0, y 2 (0) = 3. 
We choose = — 3 and f$ = —1000 in order to illus- 
trate the phenomenon of stiffness. Given the initial con- 
ditions ji(0) = 2 and 3/2(0) = 3, the exact solution is 
/3 -independent and is given by 

Exact : y\(f) = 2 exp(— t) + sin(£) , y 2 if) 
= 2 exp(-£) + cos(£). 

This example is chosen to demonstrate the performance 
of TSDM on stiff problems. We compute the solutions to 



Table 10 Results, with predictor-corrector (PreCor) and 
co = 1,for Example 3 







TSDM 




PreCor 


Simos 1998 




N 


NFEs 


CPU 


\Error\ 


CPU 


\Error\ 


NFEs 


\Error\ 


1000 


4004 


73 


1.7(-3) 


4.8 


2.9(0) 


8000 


1.4(-1) 


2000 


8004 


145 


2.5(-4) 


9.5 


4.1(0) 


16000 


3.5(-2) 


4000 


16004 


290 


2.7(-5) 


19 


3.1 (-2) 


32000 


1.1 (-3) 


8000 


32004 


584 


1.6(-6) 


38 


2.3(-2) 


64000 


8.4(-5) 


16000 


64004 


1194 


1.0(-7) 


75 


3.3(-3) 


128000 


5.5(-6) 


32000 


128004 


2546 


6.3(-9) 


150 


4.1 (-4) 


256000 


3.5(-7) 
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-•— SIMOS 



Log 10 (ErrorMax) 

-4 

-6 

-8 
-10 
-12 
-14 



5000 



Figure 4 Efficiency curves for Example 2. 



10000 15 000 



IXARU 




♦ TSDM 



NFEs 




t NFEs 

50000 100000 150000 200000 250000 

Figure 5 Efficiency curves for Example 3. 
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-\-^ — . — . — i — . — . — . — i — . — . — . — i — . — . — . — ' NFEs 
200 400 600 800 

Figure 6 Efficiency curves for Example 4. 
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- 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — CPU seconds 
5 10 15 



Figure 10 Time efficiency curve for Example 2 with predictor-corrector. 
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Log 10 (ErrorMax) 
0 





500 1000 1500 2000 2500 
Figure 12 Time efficiency curve for Example 3 with predictor-corrector. 



CPU seconds 



Example (6) with p = -3, -1000. We obtain better abso- 
lute errors than Nguyen et al (2007). This efficiency is 
achieved using fewer number of steps and less number of 
function evaluations than Nguyen et al. (2007). For exam- 
ple when p = — 3, our method generates a solution with 
error magnitude 10 -6 involving just 6 steps and 28 func- 
tion evaluations, whereas (Nguyen et al. 2007) attains the 
same error magnitude using 10 steps and 47 function eval- 
uations. When p = —1000, TSDM generates solutions 
with comparable error magnitude. We see that TSDM is 
competitive and better than the method in Nguyen et al. 
(2007) which is of order six and is thus expected to do 
better. 

An implementation in predictor-corrector mode 

In this section, we also implement our CTSDM in a 
predictor-corrector mode. The predictor is given by 

yn+i =yn + h(a 0 (u)f n ) + h 2 (k 0 (u)g„), (13) 

where 

sin(^) 



a 0 



Xq = 



2sin(^r) 



(14) 



and the corrector is given by equations (6) and (7). We 
note that when u —> 0 we use the following Taylor series 
expansion (see Simos 1998) 



0 4 

do - 1 - -g" + 120 
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+ ■ 
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A 0 — 7 loo ~T 61440 4178 
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192 ^ 61440 41287680 ^ 47563407360 



+ 



„16 



+ ... 

(15) 



As we expected, the predictor-corrector (PreCor) mode 
runs faster than the TSDM but is less accurate compared 
to the TSDM. We illustrate this by applying the predictor- 
corrector to Example 2 and Example 3. We plot the 
efficiency curves showing the accuracy versus the CPU 
computation time, and the accuracy versus the NFEs. 



Estimating the frequency 

A preliminary testing indicates that a good estimate of the 
frequency can be obtained by demanding that LIE = 0, 
and solving for the frequency. That is, solve for co given 

that (-^(^V 3) fe)+J (5) fe))) = 0, where/),; = 
2, . . . , 5 denote derivatives. We used this procedure to 
calculate co for the problem given in example (5) and 
obtained co & ±314.16, which approximately gives the 
known frequency co = 314.16. Hence, this procedure is 
interesting and will be seriously considered in our future 
research. 

We note that estimating the frequency and the choice of 
the frequency in trigonometrically-fitted methods is chal- 
lenging and has grown in interest. Existing references on 
how to estimate the frequency and on the choice of the fre- 
quency include G. Vanden Berghe et al. 2001b, and Ramos 
et al. 2010. 

Conclusion 

We have proposed a TSDM for solving oscillatory I VPs. 
The TSDM is A-stable and hence, an excellent candidate 
for solving stiff I VPs. This method has the advantages of 
being self-starting, having good accuracy with order 4, 
and requiring only two functions evaluation at each inte- 
gration step. We have presented representative numerical 
examples that are linear, non-linear, stiff and highly oscil- 
latory. These examples show that the TSDM is more 
accurate and efficient than those in Nguyen et al. 2007, 
Simos 1998, Ixaru et al. 2004, and Ozawa 2005. Details of 
the numerical results are displayed in Tables 1, 2, 3, 4, 5, 
6, 7, 8, 9 and 10 and the efficiency curves are presented in 
Figures 3, 4, 5, 6, 7, 8, 9, 10, 11 and 12. Our future research 
will incorporate a technique for accurately estimating the 
frequency as suggested in subsection "Estimating the fre- 
quency" as well as implementing the method in a variable 
step mode. 
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